############################################################################################################################################
# The Cultural Origins of Populism ----------------------------
# Yotam Margalit, Shir Raviv and Omer Solodoch  ---------------------------
# This R file contains the code necessary to replicate the analysis of Figure 5 (US case) in the main text ----------------------------
############################################################################################################################################

#### Clear environment and return memory to the operating system ---------
rm(list = ls())
gc()

#### Installing and loading needed packages -----------
# Load or install the required packages for this analysis
packages_required <- c(
  'stargazer', 'lmtest', 'haven', 'expss', 'dplyr', 'rmarkdown', 'psych', 'Hmisc','base','srvyr','broom',
  'ivmte', 'ggpubr', 'writexl', 'rio', 'patchwork', 'Lahman', 'cregg', 'xtable', 'ggmosaic','forcats',
)

# Function to check and install any missing packages
install_and_load_packages <- function(packages) {
  for (package in packages) {
    if (!require(package, character.only = TRUE)) {
      install.packages(package)
      library(package, character.only = TRUE)
    }
  }
}

# Apply the function to the list of required packages
install_and_load_packages(packages_required)


# Set Working Directory  ############
# Set the working directory where the data is stored and where to save outputs
setwd("insert here.. /JOP Replication")
plotpath<-"figures/"

#### Data ---------
data <- readRDS("data/ANES_2016_raw_data.rds")

# Trump Vote ----
data$TrumpVoter <- 0
data$TrumpVoter[data$V162058x==11]<-1
data$TrumpVoter[data$V162058x<10]<-0
data$TrumpVoter[is.na(data$V162058x)]<-NA


# Rural Resentment ----
# rural_smalltown_d
data$rural_smalltown_d <- 0
data$rural_smalltown_d[data$V167534 ==1 | 
                         data$V167534 == 2] <- 1
data$rural_smalltown_d[data$V167534==-8 |
                         data$V167534==-1  ] <- NA
data$rural_smalltown_d[is.na(data$V167534)] <- NA

# rural_resentment
data$offic_dontcare1 <- NA
data$offic_dontcare1[data$V162215 == 1|
                  data$V162215 ==2 ] <- 1
data$offic_dontcare1[data$V162215 == 3 | 
                  data$V162215 == 4 | 
                  data$V162215 == 5] <- 0


data$rural_resentment <- 0
data$rural_resentment[data$offic_dontcare1==1 & data$rural_smalltown_d==1]<-1
data$rural_resentment[is.na(data$offic_dontcare1) |
                        is.na(data$rural_smalltown_d)]<-NA
table(data$rural_resentment)

# Community disintegration ----
# Generate community_lowcontact
data$community_lowcontact<- 0
data$community_lowcontact[data$V167534!=4 &   # Non urban
                            data$V162195==2 & # community work 
                            data$V162196==2 & # meeting on school/community issue 
                            data$V162202==2 & # contacted elected local official 
                            data$V162204==2]<-1 # contacted non-elected local official 
data$community_lowcontact[is.na(data$V167534) |
                            is.na(data$V162195) |
                            is.na(data$V162196) |
                            is.na(data$V162202) ]<-NA

# Ethnocultural estrangement ---------
# White
data$race1 <- 0
data$race1[is.na(data$V161310x)]<-NA
data$race1[data$V161310x==-2|
             data$V161310x==-9]<-NA
data$race1[data$V161310x==1]<-1

# Native
data$native_2nd1 <- 0
data$native_2nd1[is.na(data$V161315)]<-NA
data$native_2nd1[data$V161315==-8|
             data$V161315==-9]<-NA
data$native_2nd1[data$V161315==1 |
                   data$V161315==2]<-1

# white_native
data$white_native <- 0
data$white_native[data$native_2nd1==1 &
                            data$race1==1] <-1 
data$white_native[is.na(data$native_2nd1)]<-NA
data$white_native[is.na(data$race1)]<-NA


 
# img_undermined_cult
data$img_undermined_cult <- 0
data$img_undermined_cult[data$V162269==1 |
                           data$V162269==2]<-1 #  somewhat / strongly agree
data$img_undermined_cult[data$V162269==5 |
                           data$V162269==4 | data$V162269==3]<-0 
data$img_undermined_cult[is.na(data$V162269)]<-NA
data$img_undermined_cult[data$V162269==-9| data$V162269==-8|  
                           data$V162269==-7|  data$V162269==-6]<- NA

# ethnocult_alienation
data$ethnocult_alienation <- 0
data$ethnocult_alienation[data$img_undermined_cult==1 &
                           data$white_native==1] <-1 
data$ethnocult_alienation[is.na(data$img_undermined_cult)]<-NA
data$ethnocult_alienation[is.na(data$white_native)]<-NA



# intergenerational backlash ---------
data$age <- data$V161267
data$age[data$V161267==-9 |
           data$V161267==-8 ]<-NA
data$age[data$V161267==90 ]<-90

table(data$age)

# Recoding birthyear into cohort
data$age_55<-0
data$age_55[data$age>54 & data$age<120]<-1
data$age_55[is.na(data$age)]<-NA
table(data$age_55)

# child_respect 
data$V162239
data$child_respect<-0
data$child_respect[data$V162239==1]<-1
data$child_respect[data$V162239==3]<-2
data$child_respect[data$V162239==2]<-3
data$child_respect[is.na(data$V162239)]<-NA
data$child_respect[data$V162239<0]<-NA
table(data$child_respect)


# child_manners

data$child_manners<-0
data$child_manners[data$V162240==1]<-1
data$child_manners[data$V162240==3]<-2
data$child_manners[data$V162240==2]<-3
data$child_manners[is.na(data$V162240)]<-NA
data$child_manners[data$V162240<0]<-NA
table(data$child_manners)

# child_obedience
data$child_obedience<-0
data$child_obedience[data$V162241==1]<-3
data$child_obedience[data$V162241==3]<-2
data$child_obedience[data$V162241==2]<-1
data$child_obedience[is.na(data$V162241)]<-NA
data$child_obedience[data$V162241<0]<-NA
table(data$child_obedience)


# child_wellbehaved
data$child_wellbehaved<-0
data$child_wellbehaved[data$V162242==1]<-1
data$child_wellbehaved[data$V162242==3]<-2
data$child_wellbehaved[data$V162242==2]<-3
data$child_wellbehaved[is.na(data$V162242)]<-NA
data$child_wellbehaved[data$V162242<0]<-NA
table(data$child_wellbehaved)


# Calculate authoritarian_index
data$authoritarian_index <- rowMeans(data[c("child_respect", "child_manners", "child_obedience", "child_wellbehaved")], na.rm = TRUE)

# Recode authoritarian
data$authoritarian <- ifelse(data$authoritarian_index >= 2.5 & data$authoritarian_index <= 3, 1,
                             ifelse(data$authoritarian_index >= 1 & data$authoritarian_index < 2.5, 0, NA))
table(data$authoritarian )

# Create older_authoritarian
data$older_authoritarian <- 0
data$older_authoritarian[data$authoritarian==1 &
                           data$age_55==1 ]<-1
data$older_authoritarian[is.na(data$authoritarian)]<-NA
data$older_authoritarian[is.na(data$age_55)]<-NA
table(data$older_authoritarian)

# Social status ---------
data$midlow_social_class <- 0
data$midlow_social_class[ is.na(data$V162128)]<-NA
data$midlow_social_class[ data$V162128==2]<-NA
data$midlow_social_class[ data$V162129==2]<-1 # Working  class
data$midlow_social_class[ is.na(data$V162129)]<-NA
data$midlow_social_class[ data$V162132==2]<-1 # Working  class
data$midlow_social_class[ is.na(data$V162132)]<-NA
data$midlow_social_class[ data$V162133==1]<-1 # Lower middle class
data$midlow_social_class[ is.na(data$V162133)]<-NA
data$midlow_social_class[ data$V162133==1]<-1 # Lower middle class
data$midlow_social_class[ is.na(data$V162133)]<-NA
table(data$midlow_social_class)

data$social_status_anx <- 0
data$social_status_anx[ is.na(data$V161342)]<-NA
data$social_status_anx[ is.na(data$V161310x)]<-NA
data$social_status_anx[data$V161342<0]<-NA
data$social_status_anx[data$V161310x<0]<-NA

data$social_status_anx[ is.na(data$midlow_social_class)]<-NA
data$social_status_anx[ data$V161342==1 &
                            data$V161310x==1 &
                            data$midlow_social_class==1]<-1 # 1. Male &   1. White, non-Hispanic & 1. midlow_social_class
table(data$social_status_anx)


##### Figure (5) Cultural Predictors and Support for Trump, ANES 2016 ------

data<-data%>%
  filter(!is.na(TrumpVoter))

data$TrumpVoter_fac<-factor(data$TrumpVoter, levels=c(0,1), labels=c("Non-Trump voter","Trump voter"))
calculate_mean <- function(var, data) {
  var_name <- as.name(var)
  data %>%
    dplyr::filter(!is.na(!!var_name)) %>%
    as_survey_design(weights = c(V160102)) %>%
    dplyr::group_by(TrumpVoter_fac) %>%
    dplyr::summarize(n = survey_mean(!!var_name, na.rm = FALSE, vartype = "ci"))
}
White_Native_img_undermined_cult_pop <- calculate_mean("ethnocult_alienation", data)
older55_authoritarian_pop <- calculate_mean("older_authoritarian", data)
rural_have_no_voice2_pop <- calculate_mean("rural_resentment", data)
social_status_anx_pop <- calculate_mean("social_status_anx", data)
distant_people_local__noturban_pop <- calculate_mean("community_lowcontact", data)

fig5_culturalstories <- bind_rows( older55_authoritarian_pop,
                                   social_status_anx_pop,
                                   rural_have_no_voice2_pop,
                                   distant_people_local__noturban_pop,
                                   White_Native_img_undermined_cult_pop,
                                   .id="Characteristic")


fig5_culturalstories$Characteristic<-factor(fig5_culturalstories$Characteristic, 
                                            levels=c(1,2,3,4,5),
                                            labels=c("Intergenerational Backlash",
                                                     "Social Status Anxiety",
                                                     "Rural Resentment",
                                                     "Community Disintegration",
                                                     "Ethno-Cultural Estrangement"))


fig5.1<-fig5_culturalstories %>%
   mutate(Characteristic = fct_reorder(Characteristic, desc(n))) %>%
  ggplot()+aes(x=Characteristic,y=n,label=n,fill=TrumpVoter_fac,colour=TrumpVoter_fac)+
  geom_col(alpha=0.9,width=0.5, position = position_dodge(width=0.7))+
  geom_errorbar(aes(ymin = n, 
                    ymax = n_upp), position = position_dodge(width=0.7),
                width = 0.5, lwd=0.4, color="gray40") +
  scale_fill_manual(values = c( "#97C1D1", "#A20021"))+
  scale_colour_manual(values = c( "#97C1D1", "#A20021"))+
  scale_y_continuous(labels=scales::percent) +
  labs(x = " ", y = "", subtitle = "", title="", 
       fill = "", color = "") +  coord_flip()+
  theme_minimal() +
  theme(legend.position = "bottom",  
        strip.text.x = element_text(size = 12),
        axis.text = element_text(size = 13),
        axis.text.x = element_text(size = 11),
        plot.background = element_rect(fill = "white", colour = "white"),
        plot.margin = unit(c(0,1, 0.5, 0), "cm"),
        panel.spacing = unit(1, "lines"))
fig5.1


fig5_coef <- data %>%
  do(tidy(lm(TrumpVoter ~                
               ethnocult_alienation+
               older_authoritarian+
               rural_resentment+
               social_status_anx+
               community_lowcontact,
             weights = V160102, data = .)))


fig5_coef<-fig5_coef%>% 
  filter(term %in% c("ethnocult_alienation", 
                     "older_authoritarian", 
                     "rural_resentment",
                     "community_lowcontact", 
                     "social_status_anx")) %>%
  mutate(term=case_when(
    term=="community_lowcontact"~"Community Disintegration",
    term=="older_authoritarian"~"Intergenerational Backlash",
    term=="rural_resentment"~"Rural Resentment",
    term=="ethnocult_alienation"~"Ethno-Cultural Estrangement",
    term=="social_status_anx"~"Social Status Anxiety"
  )) %>%
  mutate(term = factor(term, levels = c( 
    "Community Disintegration",
    "Intergenerational Backlash", 
    "Rural Resentment",
    "Ethno-Cultural Estrangement",
    "Social Status Anxiety"
    )))  # Set the order of the factors


fig_5.2 <- fig5_coef %>%  
  ggplot() + 
  aes(term, estimate, fill=term, color=term) + 
  geom_hline(yintercept=0, linetype=2, lwd=0.9, size=1.3, colour = "#ACBBBF", alpha=0.9) +
  geom_errorbar(stat = "identity", alpha = 0.5, 
                position = position_dodge(width = 0.7),
                aes(ymin=estimate - 1.96*std.error, ymax=estimate + 1.96*std.error),
                lwd=0.8, width = 0) +
  geom_errorbar(stat = "identity", alpha = 0.8, 
                position = position_dodge(width = 0.7),
                aes(ymin=estimate - 1.64*std.error, ymax=estimate + 1.64*std.error), 
                lwd=0.6, width=0) +
  scale_colour_manual(values = c("#96051A","#96051A","#96051A","#96051A","#96051A")) +
  scale_fill_manual(values = c("#96051A","#96051A","#96051A","#96051A","#96051A")) +
  geom_point(stat = "identity", alpha = 0.7,  
             size = 1, position = position_dodge(width = 0.7)) + 
  guides(color = FALSE, fill = FALSE) +
  coord_flip() + 
  labs(x = "", y = "", subtitle = "", shape="") +
  theme_minimal() +
  theme(legend.position = "bottom",  
        strip.text.x = element_text(size = 12),
        axis.text = element_text(size = 12),
        axis.text.x = element_text(size = 11),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.background = element_rect(fill = "white", colour = "white"),
        plot.margin = unit(c(0,1, 0.5, 0), "cm"),
        panel.spacing = unit(1, "lines"))

fig5<-fig5.1+fig_5.2
fig5
ggsave(fig5, #dpi = 600,
       filename = paste0(plotpath, "fig5.pdf"), device = "pdf",
       height = 4, width = 10)

ggsave(fig5,
       filename = paste0(plotpath, "fig5.tif"),
       device = "tiff",
       height = 4, width = 10, dpi = 600)